﻿
class caja
{
private:

public:
	double lbox[3];//Lado de la caja
	particula *part; //Array de particulas
	int npart;//Numero de particulas
	double phi;//Concentracion en volumen
	double rho;//Concentracion en masa
	double temp;//Temperatura de caja
	double temp0;//Temperatura de referencia
	double vol;//Volumen
	double rcut;//Distancia de corte de interaccion
	double dr[3];//Distancia a pares entre particulas
	double presion;

	///Parametros lara la link-cell (vecinos)
	int ncell;//Numero de celdas 
	int ncell2;//ncell^2
	int ncell3;//ncell^3
	double lgrid;//Lado de celda
	/*head[i] contiene el primer indice de particula contenido en celda i
	head[i]=-1 si no hay particulas en celda i*/
	int *head;
	//Contiene la lista de particulas enlazadas
	int *link;
	/*Ejemplo:
	celda i
	p1=head[i]
		- p1 es el indice de la primera particula en celda i. 
		- si no hay particulas en celda i, p1=-1
	p2 = link[p1]
		- p2 es el indice de la segunda particula en celda i
		- p2=-1 si no hay mas particulas en la celda
	etc...
	Por lo tanto, en la celda i habrá
	[p1,p2,..., pk] particulas	
	*/
	//Energias
	double epoten;//Energia potencia por particula en la caja
	double virial;
	double ekin;//Energia cinetica por particula en la caja
	//Control de la velocidad (termostato)
	double alpha; //Coeficiente de transformacion de velocidades (sqrt(Tdeseada/Treal))
	int nacc;///Paso de la simulacion
	double vel[3];///Velocidad media de las particulas de la caja

	caja()
	{

	}

	caja(int np, double xT, double xphi, double xrcut, double Dt)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		npart=np;
		temp0 =xT;
		phi=xphi;
		rcut = xrcut;
		rho = phi*6.0/pi;
		double bx = pow(pi*npart/(6.*phi),1./3.);

		for(int i=0;i<3;i++)
			lbox[i]= bx;
		part = new particula[npart];
		link = new int[npart];
		SetMeshPart();
		initForces();
		initVelocities(Dt);
		vol=lbox[0]*lbox[1]*lbox[2];

		ncell=int(lbox[0]/rcut);
		ncell2=ncell*ncell;
		ncell3=ncell*ncell*ncell;
		lgrid= lbox[0]/double(ncell);
		cout<<"Numero de celdas: "<<ncell<<endl;
		cout<<"Lado de celda: "<<lgrid<<endl;
		head = new int[ncell*ncell*ncell];
		alpha=1.0;
	}

	///Reinicia la simulacion desde un archivo
	caja(char filename[50], double Dt)
	{

		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ifstream MiArchivo (filename,ios::in);
		if (MiArchivo.is_open())
		{
			///Cargamos las caracteristicas de la caja
			MiArchivo>>npart>>phi>>temp>>rcut;
			rho = phi*6.0/3.1416;
			double bx = pow(3.1416*npart/(6.*phi),1./3.);
			for(int i=0;i<3;i++)
				lbox[i]= bx;
			part = new particula[npart];
			link = new int[npart];
			SetMeshPart();
			initForces();
			initVelocities(Dt);
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];
			alpha=1.0;
			///Cargamos las posiciones de las particulas
			for(int i=0;i<npart;i++)
				for(int j=0;j<3;j++)
				{
					MiArchivo>>part[i].pos[j];
				}
				///Cargamos las velocidades
				for(int i=0;i<npart;i++)
					for(int j=0;j<3;j++)
					{
						MiArchivo>>part[i].vel[j];
						///Construimos las posiciones anteriores de las particulas
						part[i].posOld[j]=part[i].pos[j]-part[i].vel[j]*Dt;
					}
					MiArchivo.close();
		}

	}

	~caja()
	{
		delete [] part;
	}

	///Imprime la configuracion
	void printConfig(void);
	///Coloca las particulas de forma aleatoria en una caja
	void SetRandomPart(void);
	///Coloca las particulas en red
	void SetMeshPart(void);
	///Inicializa las fuerzas
	void initForces(void);
	///Imprime la informacion de la particula considerada
	void followPart(int ind);
	///Calcula la linkCell
	void linkCell(void);
	///Calcula la lista de vecinos
	void listaVecinos(void);
	///Inicializa fuerzasK suma ewald yukawa (No probado)
	void initForcesKR(void);
	///Inicializa las velocidades y calcula posOld
	void initVelocities(double Dt);
	///Aplica condiciones contorno periodicas a particula ind
	void boundary(int ind);
	///Calcula la distancia cuadratica entre particulas
	double distancia2(int ind1,int ind2);
	///Calcula la distancia entre particulas
	double distancia(int ind1,int ind2);
	///Imprime el archivo XYZ de la configuracion de particulas
	void printXYZfile(char namefile[30], int i);
};

void caja::printConfig(void)
{
	string posPart="";
	string velPart="";
	char aux[100];
	for(int i=0;i<npart;i++)
	{
		sprintf(aux,"%f%s%f%s%f%s",part[i].pos[0]," ",part[i].pos[1]," ",part[i].pos[2],"\n");
		posPart+=aux;
		sprintf(aux,"%f%s%f%s%f%s",part[i].vel[0]," ",part[i].vel[1]," ",part[i].vel[2],"\n");
		velPart+=aux;
	}
	sprintf(aux,"%s%d%s","conf_",nacc,".txt");
	ofstream MiArchivo (aux,ios::out);
	if (MiArchivo.is_open())
	{
		MiArchivo<<npart<<" "<<phi<<" "<<temp<<" "<<rcut<<endl;
		MiArchivo<<posPart;
		MiArchivo<<velPart;
		MiArchivo.close();
	}
}

void caja::SetRandomPart(void)
{
	double pos[3];

	for(int i=0;i< npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			pos[j]=lbox[j]*double(rand())/RAND_MAX;
		}
		part[i].setPosition(pos);
		boundary(i);
	}
}

void caja::SetMeshPart(void)
{
	int np = int(ceil(pow(npart,1./3.)));
	double lx = lbox[0]/double(np);
	int n=0;
	double pos[3];
	double drift=0.5*lx;

	pos[0]=-0.5*lbox[0]+drift;
	for(int ix=0;ix<np;ix++)
	{
		pos[1]=-0.50*lbox[0]+drift;
		for(int iy=0;iy<np;iy++)
		{
			pos[2]=-0.50*lbox[0]+drift;
			for(int iz=0;iz<np;iz++)
			{
				for(int k=0;k<3;k++)
				{
					pos[k]+=(double(rand())/RAND_MAX-0.5)*0.01;
				}
				part[n].setPosition(pos);

				boundary(n);
				pos[2]+=lx;
				n++;

			}
			pos[1]+=lx;
		}
		pos[0]+=lx;
	}
	cout<<"lbox: "<<lbox[0]<<endl;
	cout<<"lx: "<<lx<<endl;
	cout<<"np: "<<np<<"  "<<int(ceil(pow(npart,1./3.)))<<endl;
}

void caja::initForces(void)
{
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].force[j]=0;
		}
	}
}

void caja::followPart(int ind)
{
	cout<<"pos: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].pos[i]<<"  ";
	cout<<endl;
	cout<<"vel: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].vel[i]<<"  ";
	cout<<endl;
	cout<<"force: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].force[i]<<"  ";
	cout<<endl;
}

void caja::linkCell(void)
{
	int mc[3], c;
	//Reset the headers, head
	for (int i=0;i<ncell3;i++)
	{
		head[i] = -1;
	}

	///Scan atoms to construct headers, head, & linked lists, lscl
	for (int i=0;i<npart;i++)
	{
		///Vector cell index to which this atom belongs
		for (int j=0;j<3;j++)
		{
			mc[j] = int((part[i].pos[j]+0.5*lbox[0])/lgrid);
			if(mc[j]>=ncell)
			{
				cout<<"mal: "<<lbox[0]<<"  "<<0.5*lbox[0]<<"  "<<part[i].pos[j]
				<<"  "<<(part[i].pos[j]+0.5*lbox[0])<<"  "<<
					(part[i].pos[j]+0.5*lbox[0])/rcut<<"  "<<ncell<<endl;
			}
		}
		///Translate the vector cell index, mc, to a scalar cell index"""
		//
		c = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
		//if(c>511)
		//cout<<"c: "<<i<<"  "<<mc[0]<<"  "<<mc[1]<<"  "<<mc[2]<<"  "<<ncell<<endl;

		///Link to the previous occupant (or EMPTY if you're the 1st)
		link[i] = head[c];
		///The last one goes to the header"""
		head[c] = i;
	}

	int ntotal=0;
	for(int i=0;i<ncell3;i++)
	{
		int cual = head[i];
		int n=0;
		while(cual!=-1)
		{
			cual=link[cual];
			n++;
		}
		//cout<<"En la celda "<<i<<" hay "<<n<<" particulas"<<endl;
		ntotal+=n;
	}
	cout<<"Numero total de particulas colocadas en Mesh "<<ntotal<<endl;
}

void caja::listaVecinos(void)
{
	double r2, rc2;
	double shield = 0.5;
	rc2=(rcut+shield)*(rcut+shield);
	for (int i=0;i<npart;i++)
	{
		part[i].nvecinos=0;
	}

	for (int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			r2 = distancia(i,j);
			if(r2<rc2)
			{
				part[i].lstVecinos[part[i].nvecinos]=j;
				part[i].nvecinos++;
				part[j].lstVecinos[part[j].nvecinos]=i;
				part[j].nvecinos++;


			}
		}
	}
}

void caja::initForcesKR(void)
{
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].forceK[j]=0;
			part[i].forceR[j]=0;
		}
	}
}

void caja::initVelocities(double Dt)
{
	//double v0=sqrt(3.0*kB*temp);
	double v02=0.0;
	double dv[3];

	dv[0]=dv[1]=dv[2]=0.0;

	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			//part[i].vel[j]=double(2.0*double(rand())/RAND_MAX-1.0)*v0;
			part[i].vel[j]=double(double(rand())/RAND_MAX-0.5);
			dv[j]+=part[i].vel[j];
		}
	}

	for(int i=0;i<3;i++)
		{dv[i]/=npart;}

	///Eliminamos el impulso inicial del sistema
	for (int i=0;i<npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			part[i].vel[j]-=dv[j];
			v02+=part[i].vel[j]*part[i].vel[j];
		}
	}

	///Transformamos las velocidades para tener la temperatura de partida deseada
	//double T = v02/(3.0*kB*npart);
	double T = v02/(3.0*npart);
	double factor=sqrt(temp0/T);
	for(int i=0;i<npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			///Reescalamos las velocidades para tener la temperatura deseada
			part[i].vel[j]*=factor;
			///Inicializamos la posicion anterior de las particulas
			part[i].posOld[j]=part[i].pos[j]-part[i].vel[j]*Dt;
		}
	}
	vel[0]=vel[1]=vel[2]=0.0;
	temp=temp0;
}

void caja::boundary(int ind)
{
	for(int i=0;i<3;i++)
	{
		int nper = int(dabs(part[ind].pos[i])/(0.5*lbox[i])); //si esta dentro de la caja 0; si esta fuera 1
		double anper = sign((double(nper)*lbox[i]),part[ind].pos[i]); //0 dentro de caja, 1 � -1 si esta fuera
		//Condiciones de contorno peri�dicas
		part[ind].pos[i]-=anper;
	}
}

double caja::distancia2(int ind1,int ind2)
{
	//Calcula la distancia entre las particulas i y j
	double dis=0;
	for (int i=0; i<3;i++)
	{
		dr[i] = part[ind1].pos[i]-part[ind2].pos[i];
		///Condiciones de contorno periodicas
		dr[i]-=int(dr[i]/(0.5*lbox[i]))*lbox[i];
		dis+=dr[i]*dr[i];
		//cout<<dr[i]<<" ";
	}
	//cout<<dis<<endl;
	return	 dis;
}

double caja::distancia(int ind1,int ind2)
{
	return	 sqrt(distancia2( ind1,ind2));
}

void caja::printXYZfile(char namefile[30], int i)
{
	char name[50];
	sprintf (name, "%s%s", namefile, ".xyz");
	ofstream MiArchivo (name,ios::out|ios::app);
	if (MiArchivo.is_open())
	{
		MiArchivo << npart<<endl;
		MiArchivo << lbox[0]<<"  "<<i<<endl;
		for(int i=0;i<npart;i++)
		{
			MiArchivo<<"O "<<part[i].pos[0]<<"  "<<part[i].pos[1]<<"  "<<part[i].pos[2]<<endl;
		}
		MiArchivo.close();
	}
}

